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ABSTRACT 

We perform three-dimensional (3D) compressible MHD simulations over many 
dynamical times for an extended range of sonic and Alfven Mach numbers and 
analyze the statistics of 3D density and 2D column density, which include proba- 
bility distribution functions, spectra, skewness, kurtosis, She-Leveque exponents, 
and genus. In order to establish the relation between the statistics of the ob- 
servables, i.e. column densities, and the underlying 3D statistics of density, we 
analyze the effects of cloud boundaries. We define the parameter space for 3D 
measures to be recovered from column densities. In addition, we show that for 
subsonic turbulence the spectra of density fluctuations are consistent with /c~^/^ 
in the case of a strong magnetic field and k~^/'^ in the case of a weak magnetic 
field. For supersonic turbulence we confirm the earlier findings of the shallow 
spectra of density and Kolmogorov spectra of the logarithm of density. We find 
that the intermittencies of the density and velocity are very different. 

Subject headings: ISM: structure — MHD — turbulence 



1. Introduction 



Magnet ohydrodynamic 



mation (see 



McKee fc Tan 



MHD) turbulence is a key eleme nt for understanding star for- 



2OO2I : iMac Low fc Klessen 1 12004| . and references therein) . In 



a new emerging paradigm, density perturbations produced by compressible MHD turbu- 
lence create molecular clouds and clumps within the clouds. This makes studies of density 
statistics very timely. In addition, recent advances in the understanding o f the properties 



of compressible magnetic turbulence (see review by ICho fc Lazarian I l2005l . and references 



therein) call for observational testing. While veloci ty statistics are currently emerging as an 
essential tool for turbulence studies (see review by iLazarian II2004I . and references therein), 
density statistics are more readily available from observations. 
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Variations of den sity have been extensively studied for a long time (see lKaplan &: Pikelner 



1970l : iDickman Ill985l ). Electron density fluctuations have become widely used to g ain insight 
into small-scale properties of interstellar turbulence (see [Armstrong et al. Ill995l . and refer- 
ence therein). Fluctuations of column density, at the same time, became a useful source o f 
information about turbulence at larger scales (IFalgarone et al. II2005I : iPadoan et al. II2003I ). 
Those studies of column densities are not lir nited to studies of spectra. Higher order struc ture 
functions have been er nployed recently fsee IFalgarone et al. II2005I : iPadoan et al. II2003I ). In 
fact, it was claimed in IPadoan et al. I (120041 ). that studies of higher order correlations could 
determine the Mach number of the interstellar flows. This motivates further studies of this 
type of statistics. 

In addition, the topology of density in the interstellar medium (ISM) can provide an 
important insight into ISM physics. For instance, the topology of the ISM is expected 
to be different for models of the ISM in w hich most of the gas is in the hot phase with 
islands of other phases embedded within (see lMcKee fc Ostriker 1119771) and in models wh ere 
a substantial part o f the volume is fi l led w ith cold and warm gas (see lCox fc Smith" 1974 ). It 



Gott et al. 



was pointed out in iLazarian et al. I f|2002r). that the topology of density can also be studied 



19861 ). However, no detailed study has been performed 



using genus statistics (see 
as far as we know. 

The issues above call for a comprehensive study of density statistics, which we undertake 
in this paper. We use an extended set of MHD numerical simulations, obtained for different 
sonic Ais and Alfvenic A4a Mach numbers, to study spectra, higher order statistics, and 
the genus of the density. We follow turbulence for many dynamical time steps (e.g. 25 
for 128^ simulations) to provide good sampling. In all cases we provide results of synthetic 
observations to compare the statistics of density with the statistics of the column density 
that are available through observations. 

In this article we investigate the scaling properties of the structure functions of density 
and the logarithm of density for compressible MHD turbulence with different sonic and 
Alfvenic Mach numbers. In ^ we describe the numerical models of compressible MHD 
turbulence. In [12] we discuss one-point probability distribution functions (PDFs). In §l]we 
study spectra of density and the logarithm of density and the anisotropy of density structures, 
and we show their dependence on Mach numbers. In ^ we perform tests of the property of 
self-similarity in the case of density structures. In ^we present the intermittency of density 
and the logarithm of density and compare the intermittency of density and column density. 
In §7] we analyze the statistical moments, skewness and kurtosis, as possible measures of 
the Ais of turbulence. In ^we study the effects of cloud boundaries, resolution and time 
averaging. In ^we analyze density structure using genus statistics. We compare the genus 
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lines to those of projected density. In ^T0\ we discuss our results and their relation to the 
previous studies. We draw our conclusions in ^ ITTl 



Numerical simulations 



We use a third-o rder-accurate hybrid essentially non oscillatory (ENO) scheme (see 
Cho fc Lazarian 1120021 . CL02) to solve the ideal MHD equations in a periodic box, 
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with zero-divergence condition V ■ B = 0, and an isothermal equation of state p = c^p, where 
p is density, v is velocity, B is magnetic field, p is the gas pressure, and Cs is the isothermal 
speed of sound. On the right-hand side, the source term f is a random large-scale driving 
force. The rms velocity SV is maintained to be approximately unity, so that v can be viewed 
as the velocity measured in units of the rms velocity of the system and B/ (Anp)^^'^ as the 
Alfven velocity in the same units. The time t is in units of the large eddy turnover time 
(~ L/6V) and the length in units of L, the scale of the energy injection. The magnetic field 
consists of the uniform background field and a fluctuating field: B = Bext + b. Initially 
b = 0. 

We drive turbulence solenoidally in Fourier space to minimize the influence of the forc- 



ated with the slow and fast modes (see 



Lithwick & Goldreich 


2001; 


Cho & Lazarian 


2003 



1 /2 

Beresnyak et al. 1 120051 . BLC05). We use units in which Va = -Bcxt/ (4vrp) ' = 1 and po = 1. 
The values of i?ext have been chosen to be similar to those observed in the ISM turbulence. 
The average rms velocity SV in a statistically stationary state is around 1. 

1 /2 

For our calculations, similar to our earlier studies (BLC05), we assume that -Bext/ (47rp) 

1 /2 

6B/{47Tpy ^6V. In this case, the sound speed is the controlling parameter, and basically 
two regimes can exist: supersonic and subsonic. Note that within our model, supersonic 
means low (3, i.e. the magnetic pressure dominates, and subsonic means high (3, i.e. the gas 
pressure dominates. 

We present 3D numerical experiments of compressible (MHD) turbulence for a broad 
range of Mach numbers (0.2 < M.^ < 7.1 and Ai^ ~ 0.7 or ~ 7.3; see Table [T]). The 
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model name contains two letters: "P" and "B" followed by a number. The letters B and P 
mean the external magnetic field and the initial gas pressure, respectively, and the numbers 
designate the value of the corresponding quantity. For example, a name "B.lP.Ol" points 
to an experiment with i?ext = 0.1 and P = 0.01. We understand the Mach number to be 
defined as the mean value of the ratio of the absolute value of the local velocity V to the 
local value of the characteristic speed Cg or Va (for the sonic and Alfvenic Mach number, 
respectively). We analyzed models with three different resolutions: 128^, 256^, and 512^ on 
a Cartesian grid with periodic boundary conditions. We drove the turbulence solenoidally at 
wave scale k equal to about 2.5 (2.5 times smaller than the size of the box). This scale defines 
the injection scale in our models. We do not set the viscosity and diffusion explicitly in our 
models. The scale at which the dissipation starts to act is defined by the numerical diffusivity 
of the scheme. The ENO-type schemes are consid ered to be relatively low diffusion ones (see 



Liu fc Osher Ill998l : iLevy. Puppo &: Russo Ill999l . e.g.). The numerical diffusion depends not 



only on the adopted numerical scheme but also on the "smoothness" of the solution, so it 
changes locally in the system. In addition, it is also a time-varying quantity. All these 
problems make its estimation very difficult and incomparable between different applications. 
However, the dissipation scales can be estimated approximately from the velocity spectra. 
In the case of our models we estimated the dissipation scale ki, at 15, 22, and 30 for the low, 
medium, and high resolution, respectively. 

We divided our models into two groups corresponding to sub-Alfvenic (-Bext = 1-0) and 
super-Alfvenic (-Bext = 0.1) turbulence. For each group we computed several models with 
different values of pressure (see Tabled]). 



3. Probability Distribution Function 

Observations usually provide 2D maps of column density in the direction of the observer 
along the line of sight (LOS). Thus, it is very important in computational astrophysics to 
study the relation between the actual 3D properties of the simulated turbulence and the 
synthetic observations prepared on the basis of these simulations. 

First, we compare the probability distribution functions (PDFs) of density and column 
density. PDFs give us information about the fraction of the total volume occupied by a given 
value of a measured quantity. For the case of compressible turbulence in which evolution is 
described by the Navier-Stokes equation and the isothermal equati on of state, the PDF of 



density obeys a lognormal form (jPassot &: Vazquez-Semadeni Ill998l ). If we select a point in 



space and assume that the density at this point results from subsequent events perturbing 
the previous density, then the final density is a product ponj(l + Si), where po is the initial 
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density at the selected point an d Sj is a small cona pression/rarefaction factor. By the power 
of the Central Limit Theorem (IKallenberg 1 119971 ). the logarithm of the resulting density, 
log(p/po) = J2i log(l + ^i)^ should obey a Gaussian distribution. 

The process of perturbation of the column density S is similar to the perturbations 
of density described above, but here a perturbation occurs only in the fraction of the LOS 
smaller than unity; thus, each individual event induces a perturbation of column density 
with a smaller r elative arnplitud e. This implies a smaller dispersion of the column density 
distribution (see lOstriker II2003I ). 



Vazquez-Semadeni &: Garcia I (120011 ) proposed the parameter rj for determining the form 



of the column density PDF in molecular clouds, which is the ratio of the cloud size to the 
decorrelation length of the density field. They defined the decorrelation length as the lag 
at which the density autocorrelation function (ACF) has decayed to its 10% level. If the 
density perturbation events are uncorrelated for rj >1, large values of the ratio t] imply that 
the Central Limit Theorem can be applied to those events. 

The decorrelation length estimated from the ACF of the density for our medium- 
resolution models ranges from about 20 cell sizes of the computational mesh for supersonic 
models to about 50 cell sizes for subsonic models, which corresponds to 5 < < 13 if we take 
the size of the computational box as the column length (in the case of medium resolution, 
it is equal to 256 cell sizes). Found values of rj signify that in our models at least partial 
convergence to a Gaussian PDF should occur. 

In the top row of Figure [1] we show PDFs of density normalized by its mean value for 
all models with medium resolution (256^). The plot on the left shows results obtained from 
sub-Alfvenic experiments {Ma ~ 0.7), and the plot on the right results from super-Alfvenic 
models (tVIa ~ 7.3). The PDF is additionally divided by the total number of mesh cells, 
which allows us to compare models with different resolutions. In the bottom row of Figure 
[T]we also see PDFs for column densities for the same models. The widths of PDF lines for 
log(S/S) are smaller then the widths of the corresponding lines for log(p/p). This confirms 
that a contribution from density perturbations to the perturbations of column density comes 
only from a fraction of the distance along the LOS giving smaller variations of the resulting 
column density. 

Our plots confirm the strong dependence of PDFs of density on the sonic Mach number, 
an alread y known and well-studied property o f density fluctu ations in compressible turbu- 



lence (see IVazquez-Semadeni fc Garcia I l200ll : lOstriker 1 12003| . and references therein) . For 



most of the models PDFs are lognormal functions. For super-Alfvenic turbulence PDFs are 
very symmetric about a vertical line crossing their maxima. However, for models with a 
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strong external magnetic field (sub-Alfvenic turbulence, AiA ~ 0.7) and very low pressure 
(supersonic turbulence, Aig ^ 1-0), the shape of the density PDFs is significantly deformed, 
and its lower value arm ends in higher densities than in the case with a weak magnetic field 
(compare models for Aig ~ 7 drawn with solid lines in the left and right plots of Figured]). 
This supports a hypothesis that the gauge symmetry for log p that exists in Navier-Stokes 
equations is broken in MHD equations because of the magnetic tension term (BLC05), which 
physically manifests itself by preventing the formation of highly underdense regions. In the 
higher density part of the distribution we do not see a similar effect, because the highly 
dense structures are created mainly due to shocks, and in this case, the magnetic tension is 
too weak to prevent the condensation. 

We show the PDFs of column density normalized to the mean value of S and divided by 
the number of pixels in the 2D map of density projected along the direction of -Bext- We see 
that the maximum and dispersion of PDFs still strongly depend on the sonic Mach number. 
Potentially this means that direct measurements of one-point statistics taken from the maps 
of molecular clouds can give us information about the speed of sound and the regime of 
turbulence taking place in these objects. See, however, discussion of the boundary effect in 

M 

In Figure [1] we included the degree of variation of PDFs in time as grey error bars. We 
note that the departure of PDFs from their mean profiles is very small in the middle part 
around the mean value. The strongest time changes are observed in the low- and high-density 
tails, but the PDFs for different models are still separable. 



4. Spectra and anisotropy of density and the logarithm of density 



The spectrum of density fluctuations is an important property of a compressible flow. 
In some cases, the spectrum of density can be derived analytically. In nearly incompressible 
motions with a relatively strong magnetic field, the spectrum of density scales si milarly to 
the p ressure as ~ if we consider the polytropic equation of state p = ap'^ (iBiskamp 



20031 ). In weakly ma gnetized MHD turbulence , however, the density spectrum is predicted 
to scale as ~ k~^^^ ( Montgomery et al. 1119871 ). In supersonic fiows, these relations might 
not be valid anymore because of shocks accumulating matter into the local and highly dense 
structures. Due to the high contrast of density, the linear relation 6p = c^dp is no longer 
valid, and the spectrum of density cannot be related to pressure so straightforwardly. In 
addition, the strong asymmetry of density fiuctuations suggests the need to analyze the 
logarithm of density instead of density itself. 
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In Figure [2] we present the spectra of fluctuations of density and the logarithm of density 
[top and middle rows, respectively) for models with different Aig- For comparison we also 
show the spectra of velocity fluctuations (Fig. bottom row). The two columns of plots 
in Figure [2] correspond to experiments with J^a ~ 0.7 {left column) and J^a ~ 7.3 {right 
column) . As expected, we note a strong growth of the amplitude of density fluctuations with 
the sonic Mach number at all scales. This behavior is observed both in sub-Alfvenic as well 
as in super-Alfvenic turbulence (see Fig. [21 top left and top right, respectively) . Although the 
spectra significantly deviate from the power law, in Table [2] we calculate the spectral index 
of density ap and the logarithm of density aiogp within the inertial ranges estimated from the 
power spectrum of velocity. The width of the inertial range is shown by the range of solid lines 
with slopes —5/3 and —2 in all spectra plots. It is estimated to be within k G (3, 20). In Table 
[2] we also show the errors of estimation which combine the error of the fitting of the spectral 
index at each time snapshot and the standard deviation of variance of Op.iogp in time. The 
slopes of the density spectra do not change significantly with A4s for subsonic experiments 
and correspond to analytical estimations (about —2.2, which is slightly less than —7/3, for 
turbulence with AiA ~ 0.7 and about —1.7, which is slightly more than —5/3, for weakly 
magnetized turbulence with Ma ~ 7.3). Such an agreement confirms the validity of the 
theoretical approximations. Nevertheless, these theoretical considerations do not cover the 
entire regime of compressible turbulence. While the fluid motions become supersonic, they 
strongly influence the density structure, making the small-scale structures more pronounced, 
which implies flattening of the spectra of density fluctuations (see values for Ai^ > 1.0 in 
Table [2]). Although, the spectral indices are clearly different for sub- and super-Alfvenic 
turbulence, their errors are relatively large (up to ±0.3). These errors have been calculated 
by summing the maximum value of the uncertainty of fits for individual spectra at each time 
snapshot and the standard deviation of time variation of spectral indices. The uncertainty 
of fit contributes the most to the total error bar (about 60%- 70%). Nevertheless, there is 
still a small probability that both spectra, for sub- and super-Alfvenic models, could overlap, 
resulting in less of a difference between these cases. 

Note that in the case of a large contrast, a more proper quantity to analyze is the 
logarithm of density; we also show spectra of the logarithm of density fluctuations (see Fig. 
[21 middle row). As we can see, the amplitude of fluctuations of logp, like the case of p, 
grows with TVI^ at all scales. However, the logarithmic operation significantly "filters" the 
extreme values of the densityjl] and they do not distort the spectra anymore. The slopes 



^As it is discussed in BLC05, the use of the logarithm of density is well justified from the point of view 
of the Navier- Stokes equations. Use of the logarithm of density provides better results than filtering with a 
threshold (see ILazarian fc Beresnyak Il2005r ). 
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of spectra of logp are much less sensitive to the sonic Mach number, which is confirmed 
by small changes of the slope with the value of Ais (see respective columns in Table |2]). 
In addition, comparing sub- and super- Alfvenic models in Table [H we see that, consistent 
with the theoretical considerations, the density gets a more shallow spectrum as the external 
magnetic field decreases (or A4a increases). For subsonic turbulence, the energy of the 
density fluctuations (and the energy of fluctuations of the logarithm of density) scales as 
~ /c"''/'^ in the case of strong Sext and as ~ k~^^^ in the case of weak -Bcxt- 

For comparison we also show the power spectra of velocity (Figure [2], bottom row). 
While the spectra of density and the logarithm of density change both their amplitude and 
slope with the sonic Mach number, for velocity fluctuations we see almost the same-shaped 
spectra within the inertial range, insensitive to the strength of sonic motions. The only visible 
differences are present in the dissipation range, where the motions are strongly influenced 
by small-scale dissipation. This dissipation depends on the pressure of the ambient medium, 
and thus on A^g, as this is the case for shocks. 

Another question is the anisotropy of density and the logarithm of density structures. In 
Figure [3] we show lines that mark the corresponding separation lengths for the second-order 
structure functions parallel and perpendicular to the local mean magnetic fieldEI In the case 
of sub-Alfvenic turbulence, the degree of anisotropy for density is very difficult to estimate 
due to the high disperion of points. However, rough estimates suggest more isotropic density 
structures, because the points extend along the line l\\ ~ For models with A4a ~ 7.3 
the points in Figure B] have lower dispersion, and the anisotropy is more like the type from 



Goldreich fc Sidhar I (jl995l . hereafter GS95), i.e. /y ~ f/^. In both Alfvenic regimes, the 
anisotropy of density does not change significantly with Ais- Plots for the logarithm of 
density show more smooth relations between parallel and perpendicular structure functions. 
The dispersion of points is very small. Moreover, we note the change of anisotropy with 
the scale. Lower values of structure functions correspond to lower values of the separation 
length (small-scale structures), so we might note that the logarithm of density structures are 
more isotropic than the GS95 model at small scales, but the anisotropy grows a bit larger 
than the GS95 prediction at larger scales. This difference is somewhat larger in the case of 
models with stronger external magnetic field (compare plots in the left and right columns 
of Figure [3]), which may signify their dependence on the strength of -Bext- The anisotropy of 
logp structures is marginally dependent on the sonic Mach number, similar to the density 



^Tlie local mean magnetic field was computed using the procedure of smoothing by a 3D Gaussian profile 
with the width equal to the separation length. Because the volume of smoothing grows with the separation 
length /, the direction of the local mean m agnetic field migh t chan ge wi th I at an arbitra. ry point. This is 
an extension of the procedures employed in ICho fc Vishi^ (|2000l) and lCho et al~l (l2002h . 
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structures. All these observations allow us to confirm the previous studies (see BLC05) that 
suggest that the anisotropy depends not only on the scale but on A^^. In comparison to the 
discussed anisotropy of density and the logarithm of density, we also show the anisotropy 
of velocity (see Fig. [31 bottom row). The anisotropy of velocity is independent of ^As too, 
but it slightly changes with Ma (compare the bottom left and right plots in Figure [3]), 
especially at small scales. For both the sub-Alfvenic and super-Alfvenic turbulence we see 
the correspondence of the anisotropy to GS95 model. From these observations we see that 
the anisotropy of the logarithm of density structures is very similar to the anisotropy of 
velocity. 



5. Extended self-similarity for density and the logarithm of density 

The statistical distribution of velocity eddies of size / is described by the moments of 
the velocity increments called structure functions 

s^^\i) = {\svm, (4) 

where 6v is a longitudinal increment 

^t;||(0 = [v(x + l)-v(x)]-/, (5) 

and (...) denotes an ensemble averaging. Note that the structure function is defined here 
with the absolute value of the argument. This ensures the positivity of the structure functions 
also for the odd orders. 

An interesting and not yet fully understood property of these functions is their scaling 
behavior. A structure function of an arbitrary order scales with the separation length I 
within the intertial range as S^^\l) ~ From such a scaling relation, one can deduce 

a self-similarity between structure functions of two arbitrary orders. Tak ing two differen t 



orders p and p', we need only a few transformations to derive a relation (see lBiskamp II2003I ) 



S(P)(1) ~ [5(P')(/)]C{P)/C(P'). (6) 

Although simple, this relation implies two important properties of the function S^'''\l). First, 
it shows how self-similarity connects the structure functions, and second, it lets us introduce 
a normalization using a structure function of an arbitrary order p'. In Kolmogorov phe- 
nomenology, the relation S^p\1) ~ Z^*-^-* has been derived analytically, implying ({p) = p/3; 
thus, the third moment has been taken as the normalization. Let us note again that the 
above considerations concern only the part of the structure function within the inertial range. 
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Benzi et al. I ( 19931 ) reported that for hydro dynamic turbulence, functions S^^\S^^^) 
exhibit a much broader power-law range compared to S^p\1). While for the inertial range 
a similarity in the scaling of two functions stems from the Kolmogorov scaling S^^\l) ~ /, 
the power-law scaling of S^p\S^^^) protrudes well beyond the inertial, into the dissipation 
gejfl This observation shows that the dissipation "spoils" different orders of the structure 



ran; 



function in the same manner. Therefore, there is no particular need to use the third moment, 
but one can use any other i noment S^^'Hl) ^ l'^^^''^ and obtain a good power law of the form 
analogous to Eq. (jS]) (see iBiskamp 1 12003| ). Though not understood, the above property 



of S^P^ allows insight into the scaling of realistically high Re fl ows with limited Re number 
simulations or laboratory experiments. In this way, the work of iBenzi et al. I (Il993l ) allowed 
much progress in studies of the properties of self-similar flows. 

All the above considerations originally concerned the properties of velocity and magnetic 
field structure. Relations between the fluctuations of density and the fluctuations of the 
other quantities like pressure, velocity and magnetic fields in the nonlinear regime of MHD 
turbulence are nontrivial to define; therefore, the scaling function and its self-similarity are 
still not well understood. Although, there is no clear analytical reason for a self-similarity 
to be valid for density structures, we check this property by calculating analogous structure 
functions 

5(^H0 = (l/(r + l)-/(r)n, (7) 

where /(r) might be an arbitrary scalar function of position, in particular, density or the 
logarithm of density. Note that we do not assume any property of scaling so far. 

We calculate structure functions for density and the logarithm of density according 
to Eq. ([7]). We point out that in general, for a normalization function one can take an 
arbitrary order of the structure function, not necessarily that resulting from the Kolmogorov 
law S^^\l) = — |e/ (where the stru cture function is defined without the absolute value 



before averaging; see iBiskamp I l2003l ) . The latter is only valid in the case of incompressible 
hydrodynamic turbulence for the velocity field within the inertial range where the viscous 
term is negligible. The difference in definition of the structure function requires additional 
comment. The averaging in both definitions is taken over the ensemble. It means that if we 
get the negative value of the increment of the quantity 6f = /(r + 1) — /(r), we can switch 
positions and take the increment 5/ = /(r) — /(r + 1) into the averaging. This operation is 
permitted, because the Kolmogorov law is valid only for the absolute values of I. 



To test if a similar scaling law exists also for density and the logarithm of density in the 



•^In practical terms this means that instead of obtaining S'^^^ as a function of I, one gets 5'-^-' as a function 
of S^^\ which is nonhnear in a way that corrects for the distortions of S''P\ 
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magnetized compressible turbulence, we present Figure HI which shows third-order structure 
functions for density, the logarithm of density, and velocity as functions of the separation 
length. The solid lines show the hypothetical Kolmogorov scaling law S^^\l) = al, where a 
is a constant. In the case of density, the models with a strong magnetic field show relatively 
large dispersion of points. In addition, the line flattens very quickly (/ > 10). In this 
situation, the region where the structure function scales linearly with the separation length 
is difficult to find. The situation is somewhat better in the case of the logarithm of density. 
Although the relations are far from linear, we can estimate the region of linear scaling 
(8 < / < 20 for sub-Alfvenic models and 5 < / < 30 for super- Alfvenic models). In Figure H] 
{bottom row) we show the third-order structure functions for velocity. We can clearly see the 
best linear scaling withing the range 15 < / < 60 for subAlfvenic models and 20 < / < 60 for 
superAlfvenic models. The linear scaling regions for the logarithm of density and velocity 
do not correspond. The dependence of third-order structure functions on the sonic Mach 
number is marginal. The only dependence can be observed in the case of velocity in models 
with a weak magnetic field. However, the sonic Mach number does not change the slope of 
the structure function. 

In Figure Owe show a test of an extended self-similarity. We plot the structure functions 
of second, fourth, and sixth order as a function of the structure function of third order. 
Plots contain functions for sub-Alfvenic models only with different sonic Mach numbers for 
density and the logarithm of density (top left and right plots in Figure O respectively). 
The alignment of points along straight lines in Figure O clearly suggests the existence of 
a self-similarity both for density and the logarithm of density. In the case of density, this 
nice picture is slightly disturbed for supersonic models, but only for S^^'^ with order p > 3. 
Structures with large density contrasts significantly distort structure functions of higher 
orders. Nevertheless, the linear relation is still well seen. Applying the logarithmic operation 
on density before computing the structure function significantly improves a linear alignment 
of points and introduces some ordering in the placement of lines (see Fig. El top right) . We 
see that lines corresponding to the structure functions of the same order are placed on the 
same line. This is especially well seen for functions of second order. This draws a hypothesis 
that the "perfect" low-order structure functions of log p have the same slope independent 
of the sonic Mach number and the structure functions are only shifted along the line of a 
fixed inclination according to the strength of the sonic motions. For higher orders the strong 
asymmetries in the distribution of log p dominate, giving the reason for the fiattening of the 
relation between structure functions of different orders. 

The existence of a self-similarity of structure functions of density and the logarithm 
of density is very useful. The question of to what degree these relations are linear induces 
us to perform a fitting, which is shown in Figure [5] {bottom row). We have chosen one 
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particular model with Aig ~ 0.7 and AiA ~ 0.7. For this model we perform linear fitting in 
the log-log domain of second, fourth and sixth order structure functions of density (Fig. [5] 
bottom left) and the logarithm of density (Fig. \5\ bottom right) as a function of the third 
order structure function of the corresponding quantity. In this case, points are regularly 
spaced and fitting on the log-log axes gives very accurate values of slopes even for high-order 
structure functions (fractional errors in the order both for density structure functions and 
the logarithm of density for the lines are presented in Fig. [5]) . Such a high degree of linearity 
suggests the existence of highly self-similar structures independent of scale, also self-similar 
in some way, because at small scales the dissipation starts to play an important role, which 
we clearly see in spectra of density, velocity, or magnetic fields. The effect of dissipation is 
not seen for relations between structure functions of different orders, which suggests that 
the dissipation acts in a consistent way at different scales. Moreover, the dissipation must 
undergo a similar scaling relation, otherwise we would get different slopes of the SFs within 
the inertial and dissipative ranges. 

Going further, we could ask if the extendend self-similarity is valid for structure functions 
of different quantities. In Figure [6] we draw plots of the structure functions of the second, 
third, fourth and sixth order of density (left) and the logarithm of density (right) as a function 
of the structure function of third order of the velocity. The existence of any similarity should 
be manifest by the linear relation here also, but as we see, this relation is far from linear. 
There are at least two reasons for this. At first we may suppose that scaling relations of the 
structure functions of a particular order, which are used as a normalization (in this paper 
we use the third moment), look different for different quantities, but it might also be the 
influence of dissipation acting differently on different quantities. Whatever the reason is, the 
problem of relations, if any exist, between the structures of different quantities involved in 
turbulence requires further investigation. 



6. Intermittency of density and the logarithm of density 

An important insight into the turbulence can be achieved by studies of turbulence self- 
similarity. This property, which is also called scale invariance, implies that fluid turbulence 
can be reproduced by the magnification of some part of it. Take the famous model of 
incompressible Kolmogorov turbulence as an example. In this model, the energy is injected 
at a large scale L and forms eddies that transfer the energy to smaller and smaller scales. At 
the scales where the corresponding Reynolds number Re/ ~ Ivi/u is much larger than unity, 
the dissipation over the eddy turnover time ti ~ l/vi is negligible. As a result, the energy 
cascades to smaller and smaller scales without much dissipation, i.e. vf/ti ~ const, which 
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gives the well-known Kolmogorov power-law scaling for the eddies of scale /, namely, vi ~ l^^^. 
The cascade terminates at the dissipation scale, which provides an inertial range from L 
to LIie~^^'^, where Re is the Reynolds number corresponding to the flow at the injection 
scale L. At the dissipation scales, self-similarity is known to fail with turb ulence forming 
non-Gaussian dissipation structures as exemplified, e.g. in iBiskamp I (120031 ). Interestingly 



enough, present-day research shows that self-similarity is not exactly true even along the 
inertial range. Instead, fluctuations tend to get increasingly sparse in time and space at 
smaller scales. This property is called intermittency. Note that power-law scaling does not 
guarantee scale invariance or the absence of intermittency. 

Intermittency is an essential property of astrophysical fluids. As intermittency violates 
the self-similarity of motions, it becomes impossible to naively extrapolate the properties 
of fluids obtained computationally with relatively low resolution to the actual astrophys- 
ical situations. In terms of astrophysics, intermittency affects turbulent heating, momen- 
tum transfer, interaction with cosmic rays, radio waves, and many more essential processes. 
Physical interpretation of intermitt ency started after the w ork by Kolmogorov, but the first 
suc cessful model was pres ented by IShe fc Leveque I (1l994j ). The scaling relations suggested 
by IShe fc Leveque I (Il994j ) related C{p) to the scaling of the velocity vi ~ l^^^, the energy 



cascade rate t," 



l- 



and the co-dimension of the dissipative structures C, 
Cip) = ^(1 - x) + C (1 - (1 - x/Cy/^) . (8) 

g 

Parameter C is c onnected to the dirnensio n of the dissipative structures D through rela- 
tion C = 3 — D (jMiiller fc Biskamp Il2000l ). In hydrodynamical turbulence, according to 
Kolmogorov scaling, we have g = 3 and x = |. Vortex filaments, which are ID structures, 
correspond to C = 2 {D = 1). For MHD turbulence we also observe current sheets, which are 
2D dissipative structures, corresponding to C = 1 {D = 2). For these two types of dissipative 
structures we obtain two different scaling relations (substituting g = 3 and x = |), 



and 



Cip) 



C(P) 



(2/3)P/3] forC 



+ 1 - (l/3)P/=^ for C = l. 



(9) 



(10) 



Relation ([9]) is often called the S 
( ITOl) the Miiller-Biskam p scaling 



against the model (see iNovikov 



Miil 



le fc Leveque scaling (IShe fc Leveque Ill994j ) and equation 



er fc Biskamp Il2000l ). There are theoretical arguments 



1994 ). but so far the She & Leveque scaling is the best for 



reproducing the intermittency of incompressible hydrodynamic turbulence. 



So far there is no similar theory for density scaling in compressible MHD turbulence. 
The presence of shocks disturbs the spectra and structure functions of density. However, as 
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we discussed in §H the logarithm of density exhibits both spectrum and Goldreich-Sidhar 
scale-dependent anisotropy (see also BLC05). 

In Figure [7] we show the scaling exponents for density {top row) and the logarithm of 
density {bottom row) for models with Aij^ ~ 0.7 {left column) and TVIa ~ 7 {right column). 
In the case of density {top row) we show only the subsonic models because of broken self- 
similarity of density structure functions in supersonic models. To overcome this problem one 
can use techniques filtering the high-density structures, which deform statistics the most. 
The application of such filtering has been discussed by BLC05. Nevertheless, we do not 
investigate this problem further in this article. When possible, we use extended self-similarity 
to obtain the exponents, applying the third-order structure function as a normalization. We 
perform such fitting procedures for all available time steps for each model. In this way we 
obtained the mean values of exponents and their standard deviation. 

Comparing plots for sub- and superAlfvenic models (Fig. [TJ top left and top right, re- 
spectively) we see that in the presence of a strong external magnetic field, density structures 
are less intermittent. Although the error bars are relatively large, the difference is undis- 
puted. Moreover, the value of the parameter C seems to be between one and two for the 
sub-Alfvenic models with very low sonic Mach numbers and approximately equal to two 
when Ais ~ A4a ~ 0.7 (see Figure [71 top left). When the magnetic field becomes weak, the 
intermittency of density grows. In addition, it grows somewhat with the sonic Mach number. 
The changes are smaller then the errors however. 

In Figure [7| (bottom row) we present the scaling exponents calculated from the structure 
functions of the logarithm of density. After this transition, the scaling exponents can be more 
precisely determined for all models, even for supersonic ones. In turbulence with a strong 
magnetic field (Fig. [TJ bottom left) we see that intermittency behaves almost identically 
for subsonic models. We can restore the scaling of density even for supersonic motions. 
The value of C is between one and two for models with A4s ^ 2.2 and exactly two for the 
model with Aig ~ 7.0. The errors in the latter case are rather small. Our models with a 
weak magnetic field reveal more intermittent structure of log p as in the case of density (see 
Fig. [71 bottom right). In these models, the scaling exponents suggest the dimension C = 2. 
Although we see some small variations of exponents with the value of Ais, they are smaller 
then the error bars. 

What we can say about all the above models is that the structure of density is intermit- 
tent and the intermittency depends on the strength of the external magnetic field. One may 
note some changes in intermittency with the strength of sonic motions, but these changes are 
weak. In general, the supersonic motions generate more intermittent structures of density or 
the logarithm of density. The only possible exception is the logarithm of density, which has 
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structure more intermittent with stronger sonic motions in the presence of a weak magnetic 
field, but with a larger uncertainty at the same time. Finally, in comparison to the density, 
the intermittency decreases when we use the logarithm of density. 

In Figure [8] we show the scaling exponents for column densities for different models. Al- 
though the estimation of exponents induce substantial uncertainties, we see a small change 
in intermittency with the sonic Mach number in all plots, especially in plots showing the 
column density integrated across B^xt in the models with a strong magnetic field. To inves- 
tigate this change deeper we show how the intermittency of column density changes with 
sonic Mach number. In Figure [9] we present the scaling exponents of the structure func- 
tions of column densities S for different experiments. As a result we obtain the relation 
of the scaling exponent and the sonic Mach number. The exponents are normalized to the 
third-order one, so Figure [3 shows exponents for the structure functions of the first, second, 
fourth, fifth, and sixth orders. We average the exponents for each model over 16 equally 
spaced moments in time t G [10,25], and we calculate their standard deviation. In this way 
we obtain also a measure of uncertainty. Figure [9] shows the exponents for column densities 
integrated along and across -Bext- In all plots we see a systematic decrease of the exponents 
for orders larger than 3 and simultaneously an increase of the exponents for orders less than 
3. This behavior signifies a progressive growth of the intermittency of density structure with 
sonic Mach number. 



7. Skewness and kurtosis 

The most straightforward statistical properties of density and column density distribu- 
tions are the mean value and the variance. However, on a deeper level we are interested to 
know the degree of asymmetry of the distribution and its relation to the normal distribution. 
These two properties are measured by the skewness and kurtosis. The mean arithmetic value 
and the variance of the distribution of, for example, density p are given by 

1 ^ 1 ^ 

i=l i=l 

where N is the number of samples or points of the mesh in the case of simulation data. 
The mean value is a less important property in our studies, so we do not consider it here; 
however, it is required to calculate higher moments. Variance measures the width of the 
distribution and is always positive by definition. 

Skewness as a measure of the asymmetry is defined through the third-order statistical 
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moment, 

The skewness for a normal distribution is zero. Negative values of the skewness indicate 
distribution that is skewed in the left direction, and positive values for the skewness indicate 
the distribution that is skewed in the right direction. By skewed left, for instance, we mean 
that the left tail is longer compared to the right tail. 

Kurtosis is a measure of whether the quantity has a distribution that is peaked or flat 
in comparison to the normal distribution. That is, data sets with high kurtosis tend to have 
a distinct peak near the mean, decline rather rapidly, and have heavy tails. Data sets with 
a low kurtosis tend to have a flat top near the mean rather than a sharp peak. A uniform 
distribution would be the extreme case. Kurtosis is defined in a similar way to skewness, 
but using the fourth-order statistical moment. 

Positive kurtosis indicates a "peaked" distribution and negative kurtosis indicates a "flat" 
distribution. 

In Figure [TD] we show the variance, skewness and kurtosis {top, middle, and bottom 
rows, respectively) for densities and column densities for all low-resolution models (models 
with resolution 128^, see Table [1]). From such a set of models we obtain the dependence 
of the moments on the sonic Mach number. Error bars are determined using the standard 
deviation of the results for each model from all subsequent steps for times t G [5, 25]. 

The variance of density fluctuations has been studied extens i vely, so its relation to 



Ms is well known flNordlund fc Padoan Ill999l : lOstriker et al. II2OOII : ICho fc Lazarian 112003 



e.g.). The variance of density grows with the sonic Mach number and consequently with the 
broadening of the PDF. The same behavior is observed in the variance of column densities. 
Lines for variances in Figure [TD] look very similar for models with strong and weak magnetic 
fields. However, the comparison of the variances for column densities integrated along and 
across the external magnetic field reveals a difference in its values, a sort of split, in the case 
of the sub-Alfvenic models (see Fig. [10], top left). In the case of super- Alfvenic turbulence, 
there is no split at all. Lines perfectly overlap. The latter is natu ral, as the super- Alfveni c 



turbulence easily tangles the external magnetic field. According to ICho fc Lazarian I (l2003l ). 
the relations between the variance of density fluctuations and the sonic Mach number are 
6p/po ~ Ais when the magnetic pressure dominates and 5p/po ~ when the gas pressure 
dominates. The former case is observed when Aig ^ A^a- Indeed, in Figure [11] we see 
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that the mean standard deviation of the density fluctuations (Sp/po) scales with A4s almost 
linearly when Ais > M.a- When M.s < M.Ai the relation is much steeper. A similar 
behavior is observed in the case of super-Alfvenic turbulence (Fig. [TTl right) although the 
gas pressure dominates in all models. For models with a very small value of the sonic Mach 
number, the relation is different. The value of {Sp/po) depends much less on JUg. 

Studies by Passot fc Vazquez-Semadeni ( 20031 ) showed the dependence of density fluc- 



tuations on the Alfvenic Mach number, which is different for different modes. When Ai^ is 
small, the slow mode dominates with the relation ~ p. The fast mode starts to contribute 
more significantly in the generation of density fluctuations with the relation ~ when 
AiA is large. A split in the a^-Aig relation for column densities parallel and perpendicular 
to Sext in the sub- Alfvenic turbulence (Fig. [101 top left) can be explained by a correlation of 
clumps and filaments of the density structure with the direction of the mean magnetic field 
resulting from the action of the slow mode as a main contributor to the generation of density 
fluctuations. In the super-Alfvenic models, the variance of the column density is independent 
of the direction of the LOS (compare to Fig. [TOl top right). Plots of variances show that 
the strength of sonic motions can be directly measured from the variance of intensity in the 
maps of molecular clouds. 

More interesting information is provided by the two other statistical measures: skewness 
and kurtosis. In Figure [10] {middle row) we see that asymmetry of density and column 
density distributions grows with the value of Aig, but only if Ais ^ 0.5. When the sonic 
Mach number falls below this value, the asymmetry of the distribution starts to grow again. 
The skewness of column density varies in the same way with the minimum placed exactly 
at the same value of Ais- In all cases, the skewness has positive values, which is interpreted 
as an asymmetry of the distribution into the direction of higher densities. The skewness is 
rather insensitive to the strength of the external magnetic field. The kurtosis, presented in 
Figure [TUl ( bottom row), oscillates around zero for models with highly subsonic motions. This 
suggests an almost perfect Gaussian distribution for density and column density. It is also in 
agreement with the fact that if density is perturbed weakly, its distribution is normal. The 
growth of kurtosis with A4s and its positive values signify a distribution "peaked" more than 
the standard normal distribution. All three measures reveal strong dependence on the sonic 
Mach number. Moreover, the utility of the above measures for the column density follows 
from the similarity to the corresponding measures of 3D density. 
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8. Additional considerations 

8.1. How does the cloud boundary influence the statistics? 

To connect our studies with column densities available from observations, we transform 
our density cube to a cloudlike shape. This operation has been done by a simple multiplica- 
tion of the density by a 3D function, which is uniform within a ball of radius R and decaying 
outside R with a Gaussian profile. In the next step we calculated column densities from the 
resulting data and performed an analysis in the same way as was previously done. In this 
simple way we obtain the synthetic map of a molecular cloud. We should emphasize here 
that the results obtained from the synthetic observations of this type are mostly applicable 
to supersonic turbulence. The restriction is due to the fact that the density fluctuations 
in subsonic turbulence are small in comparison to the mean cloud profile, so the two-point 
statistics may be dominated by this large-scale distribution. Therefore, in this section we 
narrow down our studies only to the supersonic experiments. Examples of maps obtained 
after applying the above procedure are presented in Figure [T2j The change in structure is 
visible "by eye". For highly supersonic motions (Fig. [12], left column), the column density 
exposes filamentary structure, which is created by strong shocks. The overall shape of the 
cloud is difficult to determine. In the middle column of Figure [T21 for medium sonic Mach 
numbers, the overall shape of the cloud starts to appear. The maps still show filaments but 
in a much more compact form with many spongelike elements. For subsonic turbulence, the 
cloud looks like a ball of cotton. We can recognize some whirled streaks. In this case, the 
mean shape of the cloud clearly dominates. 

To answer the question concerning the infiuence of boundaries on the statistics, we start 
by presenting the PDFs of column densities (Figure [T3l) . All models have distributions very 
different from those presented in Figure [H The PDFs are very far from the symmetry with 
respect to the mean value. The shape of the low-density part is influenced mainly by the 
boundaries of the cloud. The boundary of a cloud is an interface where the density decreases 
from the mean value of the interior of the cloud to the density of the ambient medium, so this 
explains why the shape of the low-density part looks very similar in all models independent 
of the Mach number. However, if we take a look at the high-density part of the PDFs, we 
see that here the shape does depend on the sonic Mach number. Because the densities of 
the ambient medium are too small to disturb this part of the PDF, it is clear the these are 
formed only by the internal structure of the cloud. Thus, for our cloud model, the PDFs 
from the higher density part of the cloud are more reliable. 

In Figure dHwe show three statistical moments of cloudlike-shaped densities as functions 
of moments for the original full cube data. In the six plots in Figure [THwe show variances {top 
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row), skewnesses {middle row) and kurtosises {bottom row) for sub-Alfvenic {left column) and 
super- Alfvenic {right column) experiments. We plot the values of moments obtained from 
models with different sonic Mach numbers as a function of moments for the same models but 
with densities modulated by the Gaussian distribution. The boundary does not influence 
significantly the moments if the relations presented in the plots of Figure [T3] are linearly 
increasing. On the contrary, the boundary dominates if values of the moments for Pcioud or 
Scioud do not change with the moments for Pcube and Scube- In this way we can see how the 
boundary influences the variance, skewness, and kurtosis of density and column density. 

Taking a closer look into the relation between the moments of 3D density (lines plotted 
using stars) we see that the linearity of the relation grows with the order of the statistical 
moment, which means that the boundary contributes the most in variance, less in skewness, 
and in the case of kurtosis, the contribution is almost negligible (see Fig. [TU top, middle, and 
bottom rows, respectively) . It means that if we want to avoid the influence of the boundaries, 
we should use the kurtosis as a measure of the sonic Mach number in turbulence with a proper 
scaling factor. However, we are mainly interested in the moments of 2D data, which are 
available from observations. Lines corresponding to the moments of column densities (plotted 
using triangles and squares) reveal much higher dispersion and irregularity. Nevertheless, 
we are still able to determine regions where the relation between statistical moments of 
column density calculated from cloudlike-shaped density and moments of unmodified column 
density is linear (see Fig. [T4l) . The skewness of S could also be a measure of the strength 
of the external magnetic field, because its relations, seen in middle left plot of Figure [TU 
are different for densities integrated along and across -Bcxt (lines plotted using squares and 
triangles, respectively). We do not see such differences in the results of models with weak 
Bext (see fig. [TU middle right). The relations for kurtosises in the case of S seem to be least 
helpful. The points are strongly agglomerated, with only a few points, corresponding to the 
highest sonic Mach numbers, creating a straight line. As a result of these agglomerations, 
the relations can only be used in the case of highly supersonic turbulence. 

Figure [15] shows two plots for a relation between the exponents of the structure func- 
tions of column densities for the sub- and super- Alfvenic turbulence {left and right plots, 
respectively) . these plots suggest a strong dependence of intermittency on Aig. It seems 
to be even stronger than in the plots in Figure M due to the presence of a boundary. As 
previously described, in models with low values of Aig, the contribution in the structure 
functions mainly comes from the large-scale structure of the cloud. In this case, the ex- 
ponents are almost independent of Aig (see the subsonic part of the plots in Figure [T^ . 
However, when the sonic motions become stronger, they are able to produce a high con- 
trast in density structures, enough to dominate in the structure functions. Thus, we are 
able to study the intermittency from the intensity maps only in turbulence with sufficient 
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values of Aig- However, to overcome this problem in the case of subsonic turbulence, we 
should restrict analogous analysis only to the intensity maps of the interior of the cloud. In 
both cases, we should treat the statistics obtained from the intensity maps very carefully, 
including the boundary effect as much as possible. 

8.2. How do the statistics of density depend on resolution? 

One could ask how our results depend on the resolution of our calculations. Undoubt- 
edly, in models with a higher resolution the small-scale structure of density is better resolved. 
We have seen in previous sections that these small-scale structures can significantly affect 
the spectra and structure functions. In this section we examine to what extent our results 
depend on numerical resolution. 

In Figure [12] we show PDFs and spectra of density for experiments with the same 
physical conditions, but for three different resolutions: 128^, 256^, and 512^. As a reference 
example we show only models with Aig ~ 0.7 and A^a ~ 0.7. Plotted functions for the low 
and middle resolutions are the result of averaging over several time steps to get the results 
independent of time variations. For the highest resolution we have only a single snapshot in 
time. Figure [TBI (left) suggests that the distribution of density is almost insensitive to the 
resolution. We see very similar shapes and values of lines. The only visible difference is a 
tail at higher densities for the high-resolution experiment, which can be a result of time- 
dependent fluctuations of the highly compressed small-scale structure. It could also be a 
result of better resolved shocks due to the smaller numerical diffusivity in the case of the 
high-resolution model, which smooths the density compression at shocks less, resulting in a 
higher density contrast. A progressive growth in the numerical diffusivity is especially well 
seen in the spectra of density (see Fig. [16], right). Due to a high diffusivity, the inertial range 
is very short for the low-resolution model, but its slope is consistent with slopes for higher 
resolution experiments. Almost perfect covering of the dissipation ranges for all resolutions 
in the right plot of Figure [TB] shows that the numerical dissipation acts in the same way on 
the small-scale density structure independent of the resolution. The resolution changes only 
the length of the inertial range. 

Figure [17] shows the scaling exponents for density for three different resolutions. To 
estimate exponents and errors we calculate the mean value and standard deviation from 
results for all available time snapshots at t > 5 for a given experiment. Figure [T7] suggests the 
growth of the intermittency of density with resolution. The intermittency is more sensitive 
to the resolution. It is especially well seen at higher orders p. Simultaneously, we see bigger 
errors for lower resolutions. But bigger errors could also indicate that the intermittency of 
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structure is oscillating stronger over time. An interplay of forces creating small-scale and 
dense structures through shock formation and diffusion (also numerical) and destroying such 
structures could establish an equilibrium of some sort when the development of turbulence 
reaches a stationary state. In this kind of analysis, there are also statistical restrictions due 
to the sampling size effects. Structure functions are reliable to some order depending on the 
resolution of numerical models. In our studies we have the size of analyzed samples ranging 
from about 2-10^ to over 10*. In such a case, the maximum order should be for p = 8. This 
explains the increase of the error bars. Our conclusion is that we should expect a somewhat 
higher intermittency of density in experiments with higher resolution. 

In some sense our study sends a warning to the low-resolution studies of the She-Leveque 
scaling. We see that while the extended self-similarity may be fulfilled for the low-resolution 
studies, their higher order scalings can still be different from higher resolution runs. 

8.3. How do the statistics of density depends on time averaging? 

In our analysis we used 16 snapshots of data uncorrelated in time for the low-resolution 
models and from eight up to 17 snapshots for the medium-resolution models depending on 
the modelE Our results show a relatively significant magnitude for the errors (see Figures [TJ 
O and [TU]) understood here as the measure of the time variance of the analyzed quantities. 
Continuous interaction between waves results in creation of the structures. We observe 
smooth transitions from one structure to another. The process is stochastic. It means that 
as with all stochastic processes, it should not be considered for one particular snapshot, but 
a relatively large number of snapshots should be taken into account. For instance, if we 
had used only one snapshot, we would come to the wrong conclusions about whether the 
intermittency is higher or lower for a particular case or whether the statistical moments 
scale differently with A4s- This is especially important for measures with a relatively high 
level of uncertainty like scaling exponents or statistical moments for the column densities. 
In addition, in numerical simulations the time step is defined by the maximum characteristic 
speeds in the system. These speeds depend on the local amplitudes of the velocity and 
magnetic field, which can be different for different models. In this case, a small shift in time 
at which we compare snapshots for different models would result in completely different 
conclusions for measures strongly varying in time. In this way, the analysis of turbulence for 
only small time averaging produces results that could be very unreliable. 



"^The accurate number of snapshots in each model is given in Table [TJ 
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9. Topology of density structures 

To study the topology of density structures, we applied the genus statistic s that are 



successfully used in the cosmic microwave background studies ( IGott et al. Ill986l . for exam- 
ple). For a given density threshold we define an isosurface of density separating higher and 
lower density regions and evaluate the genus of this contour. Roughly speaking, the genus 
of an isosurface is the number of holes in this surface minus the number of isolated regions 
for a given threshold. This difference changes with the change of the threshold. We plot the 
genus of these contours as a function of the density threshold. 



To perform genus statistics we use a Fortran algorithm published by IWeinberg I (119881 ). 
From the analysis of the genus of isosurfaces of density, we can extract information about 
the density distribution and about the topology of structures that are created by turbulence. 
The results and plots in this section refer to the medium-resolution models. 

In Figure [TS] [left column) we present the distribution of genus for three different sonic 
Mach numbers, one for the subsonic case (Fig. [181 ^op left) and two for the supersonic cases 
(Fig. [iHl moddle left and bottom left). For the subsonic case in which the distribution is 
almost symmetrical with respect to the mean density value, the shape of the genus line sug- 
gests a spongelike density structure with occasional overdense clumps or under-dense holes. 
We observe a somewhat higher number of higher amplitude dense regions than rarefactions, 
because on the high-density side of the distribution, the line decreases to highly negative 
values (down to about -80). Correspondingly the low distribution of rarefactions reaches a 
value of only -50. For supersonic models with Aig around 2.5 or 7.0, the genus line looks 
much different. Although the structure seems to be also spongelike (the distribution in the 
vicinity of the mean density value is broad and highly positive), the low-/3 turbulence pro- 
duces an extended tail in the high-density direction. This suggests the existence of many 
highly dense singular structures. It may be evidence for the formation of strong shocks, 
which perturb and compress density locally in small volumes (see BLC05). Density in these 
regions grows to a value 25 times larger than the mean value in the case of = 7 turbu- 
lence. In connection with this we observe the disappearance of the low-density tail in the 
distribution, which means that underdense holes are short lived. 

We can ask the question about whether the genus statistics of integrated density along 
one chosen axis can give us adequate information about the 3D topology of the density 
structure. With such a relation, we could interpret maps of column densities obtained from 
observations of the diffuse ISM and molecular clouds. In Figure [18] {right column) we show 
the genus statistics of density integrated along the X-axis (or along the external magnetic 
field, which is also direc ted along the X-axis) . We applied the 2D version of the algorithm 



for genus computations (IMelott Ill990t ). A 2D genus is interpreted as the number of isolated 
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2D high-density regions minus the number of isolated 2D low-density regions. While the 
3D genus is symmetric about the mean value, the 2D version is antisymmetric. For the 
high-/5 case (Fig. [TS|, top right) the genus distribution is very symmetrical, and we only 
see one rarefied hole structure with column density lower than 180. This corresponds to 
the Gaussian distribution of the integrated density for this model. For the low-/? cases the 
structure is extended more into high-density structures, hke in the 3D case of the genus line, 
and we do not observe any symmetry there. Although the number of more dense regions is 
somewhat lower than the less dense regions, the maximum is shifted to the right. We can 
say that all the dense clumps have densities above 300, while holes have densities below 150. 

It was previously discussed that clumps with density much higher than the mean value 
can distort the power spectrum and hide the structure of the density created by motions at 
small scales (BLC05). Following BLC05 we applied these techniques to test how the genus 
statistics would be affected. In Figure [19] we present the genus lines for the logarithm of 
density {left column) and the corresponding 2D genus lines for the logarithm of the integrated 
density {right column). From the 3D genus plots (Fig. [19], left column) we see that the lines 
are almost symmetric about zero. The most non symmetry reveals itself in the supersonic 
models, where the maximum shifts to the negative values of log p. For integrated density 
plots (fig. [191 right column), we do not see any significant differences between the results 
for different Mach numbers. For high Mach numbers, the logarithm of density has a higher 
extremum, which is not unexpected. 

We have tested the genus statistics for density with the cutoff of high-density peaks. 
We do not show the plots for them, because the procedure of cutting off the density above 
some level simply does not change the line in plots of Figure [18] below the cutoff value and 
only removes the high-density tail. 

The interesting information would be the sonic Mach number retrieved from the genus 
statistics of column density maps. The sonic Mach number can be obtained from the length 
of the high-density tail of the genus line. It is just another way to measure the maximum 
value of density, which is strongly related to the value of A^^, which is explained more in ^ 
However, more important information that the genus can give is the type of structure. For 
low sonic Mach number models we have both underdense holes and overdense clumps with 
approximately the same count of each. This is supported by the top right plot in Figure [TU 
where we observe distribution a roughly symmetric about the mean value. However, with the 
increase of sonic Mach number, turbulence produces more high-density structures, extending 
the right tail of the genus line significantly. This means that the structure starts to look 
rather like a medium with lower density in the most of volume, but with many singular ID 
or 2D structures with much higher density. The reader can compare this picture to Figure 
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[T2l which supports our claims. We can conclude then that the genus statistics can be a 
valuable tool to study the structure of molecular clouds. 



10. Discussion 



10.1. Density in MHD turbulence 

We feel that it is important for studies of turbulent density to be understood in the 
context of broad MHD turbul ence. There has lori g been an understanding that MHD 
turbulence is anisotropic (e.g. IShebalin et al. 1 1198311. A ve ry important insight into in- 



compressible MHD turbulence by iGoldreich fc Sidhar I (Il995l ) (GS95) has been followed by 



Cho &: Lazarian 



progress in the understanding of compressible MHD t urbulence (CL02 

2OO3I : iLithwick fc Goldreich lEooH : IVestuto et aPboOsh . In particular, simulations in CL02 
showed that the Alfvenic cascade evolves on its own and the exchange of energy between 
Alfven, slow, and fast modes is marginaljfl This understanding differs from some earlier 
claims of a high degree of coupling of compressible and incompressible motions. Inciden- 
tally, the notion of high coupling was used t o explai n the fast damping of co mpressible MHD 
turbulence reported in IStone et al. I (119981 ) and in iMac Low et al. I (119981 ). In view of the 
current understanding of MHD turbulence we clai m that the fast decay o f MHD turbulence 



is due to the fast decay of the Alfvenic m ode (see ICho fc Lazarian 1 12003| ) , the effect that is 



present even in incompressible MHD (see ICho et al. 



2OO2I 



In view of the aforementioned picture of MHD turbulence, density perturbations arise 
from the slow and fast modes. Alfven modes, however, shear and cascade the slow modes 
(GS95) and therefore affect the density spectrum. Note that for supersonic turbulence the 
slow modes get modified and produce shocks [the spectrum E{k) ~ k""^]. The large-scale 
shocks produced compress density and produce a shallow spectrum, as we observed in the 
simulations above (see also BLC05). 

In terms of the study of PDFs and spectra, our present work goes beyond the one in 
BLC05, as we provide a wide r parameter space su rvey. However, we reach conclusions similar 
to those in BLC05 (see also lPadoan et al. II2004I ). Interestingly enough, our results are also 



The expression proposed and tested in ICho fc Lazarian 1 12003[ ) shows that the couphng of Alfvenic and 
compressible modes is appreciable on the injection scale if the injection velocity is comparable to the total 
Mach number of the turbulence, i.e. with {V^ + c^)^^^, where Va and Cs are the Alfven and sound velocities 
respectively. However, the coupling becomes marginal at smaller scales, as turbulence cascades and turbulent 
velocities get smaller. 
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similar to those obtained in non- magnetized fluids by iKim &: Ryu I (120051 ) . Indeed, we find 
that the spectrum of density gets shallower than the Kolmogorov one as the sonic Mach 
number increases. The reason for this is the formation of small high-density compressed 
regions that dominate the spectrum. If the logarithm of density is used, these regions do 
not dominate and the underlying spectrum infiuenced by anisotropic Alfven mode shearing 
emerges. As for the PDFs, they are close to lognormal for sub-Alfvenic turbulence even if 
the turbulence is supersonic. For super-Alfvenic turbulence the PDFs can be approximated 
by a lognormal distribution only for low sonic Mach number. 

The remarkable regularity of log p should have a physical explanation. A possible one is 
related to the multip licative symmetry with respect to de nsity in the equations for isothermal 
hydrodynamics (see iPassot &: Vazquez-Semadeni I Il998l ). This means that if a stochastic 
process disturbs the density, it results in perturbations of density being multiplied rather 
than added together. Consequently, the distribution of density for a Gaussian driving of 
turbulence tends to be lognormal, which is consistent with our PDF measurements. Magnetic 
forces should affect the multiplicative symmetry above. However, they do not affect the 
compressions of gas parallel to the magnetic field. Compressions in magnetically dominated 
fiuid will be of the highest intensity and therefore most important. They are sheared by 
Alfvenic modes, as their own evolution will be slower than that imposed by Alfvenic cascade 
(see theoretical considerations in GS95. ICho fc Lazarian II2003I ). We note that shearing does 
not affect the PDFs, but it does affect the spectra and the anisotropy of the turbulence. 
This explains the close relationship between the properties of velocity in our simulations and 
those of density. 



10.2. Properties of 3D density 

We confirm the existence for density of a rather mysterious property of turbulence called 
"extended self-similarity". While our spectra have a limited inertial range, the moments of 
density normalized by the third moment show nice power laws. Interestingly enough, we 
find that the extended self-similarity is not true if one uses the normalization over the third 
moment of velocity. While we still may wonder about the origin of this phenomenon, the 
extended self-similarity is an indispensable tool for studyint turbulence intermittency. 

Our study of the intermittency of density and the logarithm of density reveals the de- 
pendence of the She-Leveque exponents on both sonic and Alfven Mach numbers. Moreover, 
we note a rather weak dependence on the intermittency when < 1. For supersonic tur- 
bulence the dependency is much stronger. The intermittency changes with the strength of 
the external magnetic field may frequently be more important. We observe less intermittent 
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structures when B^xt is strong (see Figure [7]). 

The question of anisotropy of structures is in a natural way related to the presence of 
magnetic field. When no magnetic field is present, there is a lack of any force that could 
organize these structures. On the other hand, the introduction of the external magnetic field 
to the turbulent system induces the Alfvenic mode, which is able to imprint its structure 
on random densities. In MHD turbulence, one direction is chosen, the direction of B^xt, 
so the anisotropy should be observed. In our sub-Alfvenic simulations the anisotropy of 
density is difficult to distinguish due to large dispersion of the structure functions. The 
situation improves in the case of models with a weaker magnetic field. Here the results 
are in accordance with the GS95 model (i.e. l\\ ~ ^l''^)- The anisotropy of the logarithm 
of density and velocity structures show similar GS95 behavior in all models with observed 
dependence on the scale and independence of J^s- However, the structures of logp and 
velocity for sub-Alfvenic models seem to be less anisotropic at the small scales. 

Some basic properties of the density distribution in turbulence have been already stud- 
ied. However, most papers were focused mainly on the variance of distribution. We analyzed 
here also two higher order statistical moments: skewness and kurtosis. We observe a strong 
dependence of these moments on A^^. For instance, skewness grows with the sonic Mach 
number for mid- and supersonic turbulence. For highly subsonic turbulence we observe the 
growth of skewness but with a decreasing A4s- Kurtosis in this case is zero and starts to 
grow rapidly when the value of Aig exceeds about 0.7. Combining these two measures can 
allow one to determine the value of Aig in turbulence. We do not observe any significant 
dependency of these moments on A^a- 

Genus statistics reveal the changes of the structure of density with sonic Mach number. 
For lower values of Aig density has a normal distribution, signifying roughly the same num- 
ber of separate rarefactions and condensations. For higher sonic Mach numbers, overdense 
structures start to dominate over the underdense ones. This changes the structure of density 
significantly. For low sonic Mach numbers, the structure of the density is spongelike. In 
the supersonic turbulence, the shocks produce separated high-density ID and 2D structures 
embedded in the medium of much lower density. In this sense, the genus statistics confirm 
the topology with the intermittency. 



10.3. What can be found from observations? 



Many surveys provide extensive information on column density statistics. Those include 
Hi surveys, e.g. the GALFA (http : / /www . naic . edu/alf a/ galf a/ ) survey of galactic Hi 
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distribution. The main goal of another observational project, the Wisconsin H-Alpha Mapper 
(WHAM) , is to provide a survey of H a emission from the ISM over the entire northern sky 
( Haffner et al. lEo03 : Hill et aLlboOS ). The observations provide the dispersion and emission 
measures, which can be interpreted as the column density and the integral of the square of the 
electron density along the entire LOS. The survey supported by the numerical studies helps 
in better understanding the density and distribution of diffuse plasma in a major component 
of our Galaxy, the warm ionized medium. Combining data from different surveys and using 
the techniques in this paper, it is possible to get valuable insight into processes taking place 
in different ISM phases. 

If the observab le column density distributions are concerned, our study confirms the 
earlier conclusion in lOstriker I (120031 ) that the PDFs are affected by the LOS distribution of 
density, and therefore, the relation between the column density PDFs and the underlying 
3D PDFs is rather involved. On the contrary, the relation between the 3D and 2D density 



spectra is very straightforward, as was shown earlier by many authors (see IStutzki I Il999 



Lazarian Ill995l ). This makes the observable column density spectrum a good measure of the 
underlying 3D density spectrum if the value of the spectral index is concerned. However, 
the anisotropies of the density distribution that can be revealed for low sonic Mach number 
by density itself and for high Mach number by the distribution of log p are affected by the 
LOS integration that obtaining the column densities involves. The problem here is that the 
GS95 anisotropies of the magnetic field and velocities that induce anisotropies of the density 
are defined in the s ystem of reference aligned with the local magnetic field (see discussion 



m 



Cho et al. 1120031 ). LOS integration averages the local field directions, and therefore, only 



the anisotropies have to be calculated with respect to the mean magnetic field. This both 
decreases the degree of anisotropy and erases the scale dependence of the anisotropy. The 
residual anisotropy is the anisotropy of the turbulence at the injection scale. Therefore, 
for super- Alfvenic turbulence for which magnetic fields are dynamically unimportant at the 
injection scale, the anisotropy of the column density statistics vanishes altogether. 



Padoan et al. I (120031 ) studied maps of CO intensity of molecular clouds in Taurus and 



Perseus in terms of structure functions. The integrated int ensity maps can b e translated into 



column densities by simple multiplication by a factor (see lOssenkopf II2002I ). iPadoan et al 



(120031 ) show that the scaling exponents of the structure functions, normalized to the third 
order, follow the velocity scaling of the supersonic turbulence. This suggests that the den- 
sity scaling can be de scribed by a similar hierarchical model, like the one proposed by 
She fc Leveque I (119941 ). Our study shows that the intermittencies of densities and velocities 
are intrinsically different and one should exercise care in relating the relevant She-Leveque 
exponents. Although the intermittencies can reflect the sonic Mach number, this dependence 
that we observe is weak for subsonic models and the changes are below the errors (see Figure 
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|9]). The low- resolution models showed the weak dependence of the intermittency of density 
on the sonic Mach number. However, the intermittency is very sensitive to the resolution, 
which was stressed in §8.21 The growth of error bars with the order of exponent p is due 
to the amplification of any irregularities in the structure function by higher order moments. 
Another difficulty is the change of the intermittency while we switch from density to column 
density, as seen by the comparison of the scaling exponents in Figures [7] and [H We note, 
that the intermittency of the column densities in general is lower than the intermittency of 
density. In such situations we can only speak about some rough predictions taken from the 
observational data. 

In addition, we see that the effects of cloud boundaries change the intermittency itself 
depending on Aig. For instance, the highly subsonic and supersonic turbulence reveal more 
intermittent structure after including the boundary. The effect is especially strong for tur- 
bulence with high Ais and grows with its value. The strength of the external magnetic field 
does not infiuence significantly the effects of the cloud boundary. 

We also present a number of different measurements that can be helpful in the determi- 
nation of the properties of density structures in real cosmic objects, like molecular clouds. 
We compared the statistical measures applied for the 3D density and corresponding 2D col- 
umn densities. By such comparison, we were able to regain information about the spatial 
structure from the synthetic maps of the projected density. For instance, we have shown, 
that the PDFs of density and column density can be related in particular cases. 

We find that the statistical moments skewness and kurtosis, which are easy to obtain, can 
be very useful. They can be employed to determine the sonic Mach number for supersonic 
turbulence. The subsonic regime of turbulence should be treated more carefully, because 
the cloud boundary effects might be very important in this case, as was shown in §8.11 
Boundary effects influence the PDFs and the scaling of density structure (see Figures [12] 
and [151 respectively). In the case of PDFs, however, some information can still be regained, 
because the boundary of the cloud mostly influences the low-density tail of the distribution. 
Three-dimensional density spectra can be recovered successfully from observations of column 
densities. 

Another promising technique for the topology of density and column density study 
is genus statistics. Although it is a completely different technique than one- or two-point 
statistics, we have shown that it can both help us with understanding the topology of density 
structures in turbulence and relate the results obtained for density and column density (see 

Some measures presented in the article, like statistical moments and spectra, should 
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be easily obtained from the molecular maps. Such observation of quantities can be directly 
compared with our numerical results, allowing quick determination of the conditions in 
investigated maps. Other measurements, like scaling exponents and genus statistics require 
more effort to calculate, but can provide more information related to the density structure. 
Thus, the techniques of the study of turbulence using 2D maps of observed column densities 
are very useful. 

11. Summary 

In this paper we investigated the statistics of density structure in compressible MHD 
turbulence using numerical experiments. We analyzed PDFs, spectra, structure functions, 
and the intcrmittency of density and the logarithm of density for experiments with different 
sonic Ais and Alfven AiA Mach numbers. We addressed the issue of whether the numerical 
results can be compared with observations. We found: 

• The amplitude of density fluctuations strongly depends on Ms both in weakly and 
strongly magnetized turbulent plasmas. 

a) PDFs of density fluctuations have a lognormal distribution. 

b) The width of the PDF of logp grows with Ms- 

c) The flattening observed in density spectra is due to the contribution of the highly 
dense small-scale structures generated in the supersonic turbulence. 

• Fluctuations of the logarithm of density are much more regular than those of density. 
The logarithm of density exhibits the GS95 scalings and anisotropics. The PDFs are 
close to Gaussian. 

• Structure functions of density and the logarithm of density reveal self-similarity. The 
density field is not self-similar to the velocity one. 

• The intcrmittency of density structure depends on both Ms and Ma- We report only 
weak dependence on for subsonic turbulence. In the case of supersonic models the 
dependence is much stronger. 

• Additional care about resolution is necessary for intcrmittency studies. The She- 
Leveque exponents results can be affected by the resolution. 

• Similar to their 3D counterparts, the 2D measures of column densities, i.e. the variance, 
skewness and kurtosis, depend on Ms and can be used to determine the value of Ms 
from observations. 
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• The cloud boundary effects are important and should be accounted for by the study. 
The most sensitive to the effects of cloud boundaries are the PDFs and intermittency 
studies. 

• The genus statistic is sensitive to A^^. The 2D genus analysis is a promising tool for 
studying the topology of the ISM. 
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Fig. 1. — Normalized probability distribution functions (PDFs) of density p/p {top row) 
and column density S/E integrated along the line of sight parallel to B^xt {bottom row) for 
models with different values of A4s and for J^a ~ 0.7 {left column) and A4a ~ 7.3 {right 
column). Grey error bars signify the standard deviation of the PDFs, showing the strength 
of the departure at a single time snapshot from its mean profile averaged over time from 
t > 5 to the last available snapshot. In this plot we include all models listed in Table [1] 
computed with medium resolution (256'^). 
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Fig. 2. — Spectra of density {top row), the logarithm of density {middle row) and velocity 
{bottom row) for experiments with different values of M.s and with A4a ~ 0.7 {left column) 
and M.A ~ 7.3 {right column) for models with medium resolution (256^). Grey error bars 
signify the variance of spectra with time. The solid lines with slopes —5/3 and —2 cover the 
inertial range used to estimate the spectral indices of p and log p spectra. 
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Fig. 3. — Anisotropy for the second-order structure function for density {top row), the 
logarithm of density [middle row), and the velocity field [bottom row) for experiments with 
different values of Ais and with A4a ~ 0.7 [left column) and A4a ~ 7.3 [right column) for 
models with resolution 256^. 
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Fig. 4. — Scaling law for density, the logarithm of density, and velocity {top, middle, and 
bottom rows, respectively) for medium-resolution experiments with AIa ~ 0.7 {left column) 
and A4a ~ 7.3 {right column). 
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Fig. 5. — Top row: Tests for some sort of self-similarity for density {left) and the logarithm 
of density (right). Plots show structure functions of three different orders p=2, 4, and 6 as a 
function of the third-order structure function for models with Ai^. ~ 0.7. Structure functions 
of models with different Aig form lines in different ranges of the plot. The linear shape of 
the lines suggests the existence of self-similarity. For high values of Aig, the dispersion of 
points is higher, but the points are still aligned along a line. In the case of logp {top right) 
the structure functions of lower orders preserve the value of the slope for a wide range of 
sonic Mach numbers. Strong density fluctuations disturb the slope of higher order structure 
functions, bottom row: Structure functions SFp{1) of density (left) and the logarithm of 
density (right) for three different orders p=2, 4, and 6 as a function of the third-order 
structure function for models with AiA ~ 0.7 and Aig ~ 0.7. Points represent the computed 
values of structure functions, and lines are a linear fit to the points. The fitting procedure 
gives the value of the slope with very small uncertainty (less than 1.5% for density and less 
than 0.4% for the logarithm of density). 
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Fig. 6. — Structure functions of density {left) and the logarithm of density {right) as a 
function of the third-order structure function of velocity. The relation is not linear, so 
there is no simple similarity between density and velocity structure functions. An important 
consequence of this nonlinearity is an arbitrary normalization, which is not because of the 
Kolmogorov law relating third-order structure functions with the separation length I. 
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Fig. 7. — Scaling exponents for density {top row) and the logarithm of density {bottom 
row) normalized to the third-order exponent for models with M.a ~ 0.7 {left column) and 
A4a ~ 7.3 {right column). For density only the values for subsonic experiments are shown 
from Ais ~ 0.2 (triangles) to ~ 0.7 (stars) for models with medium resolution. For the 
logarithm of density the values for subsonic (triangles, diamonds, and stars) and supersonic 
(squares and crosses) experiments arc shown. Lines show theoretical scalings: dotted for 
Kolmogorov and dash-dotted and dashed for She-Leveque with parameter C = 1 and 2, 
respectively {g — S and x — ^). 
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Fig. 8. — Scaling exponents of the column density integrated along and across the direction 

of i?cxt {top and bottom rows, respectively) for experiments with A4a ~ 0.7 and ~ 7.3 {left 
and right columns, respectively) for models with medium resolution. We use the third-order 
exponent to normalize. Lines show theoretical scalings: dotted for Kolmogorov and dash- 
dotted and dashed for She-Leveque with parameter C = 1 and 2 respectively {g = 3 and 
x=l). 
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Fig. 9. — Scaling exponents of the structure functions of column density E projected along 

{top row) and across {bottom row) the external magnetic field Sext for sub-Alfvenic {left 
column) and super- Alfvenic {right column) turbulence for models with low resolution (128^). 
We show the exponents for the structure functions from the first to the sixth order excluding 
the third one. 
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Fig. 10. — Statistical measures describing the distribution of density p and column densities 
S projected along and across the external magnetic field i?cxt for models with M.a ~ 0.6 {left 
column) and M.a ~ 7.2 {right column) for models with low resolution (128'^). Plots in the 
top row show variances of the density distribution. In the middle row we show skewness and 
in the bottom row kurtosis. In each plot we combine lines for densities with a corresponding 
scale on the left axis of the graph, and lines for the column densities with the scale on the 
right axis of the graph. Solid lines signify a polynomial of fourth order fitted to the points. 



-44- 




Fig. 11. — Relation between the mean value of the standard deviation of density fluctuation 
(Sp/po) and the sonic Mach number Aig for models with low resolution. Two lines show 
the analytical relations derived for magnetic-pressure-dominated turbulence (~ Adg, dashed 
line) and the gas-pressure-dominated turbulence (~ M^, dot-dashed line). 
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Fig. 12. — Examples of generated synthetic maps for high-resolution models with the same 
strength of the external magnetic field and different sonic Mach numbers {A4s = 8.3, 2.1, 0.6 
for left, middle and, right columns, respectively) . Plots in the top row show column densities 
integrated along -Bext, and in the bottom row, densities are integrated across -Bcxt- 
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Fig. 13. — PDFs of density p shaped in the form of a cloud {top row) and with its column 
densities E projected along and across the external magnetic field B^^i {middle and bottom 
rows, respectively) for sub-Alfvenic and super- Alfvenic turbulence {left and right columns, 
respectively). We used models with medium resolution here. 
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Fig. 14. — Comparison of the variance, skewness, and kurtosis {top, middle and bottom rows, 
respectively) of the cloudhke-shaped density plotted as functions of the unshaped full cube 
densities. In the left column we present results for sub-Alfvenic turbulence, while the right 
column contains results for super-Alfvenic turbulence. Each plot consists of the moments 
for models with different sonic Mach numbers for density and column densities integrated 
along and across the external magnetic field. The low-resolution models are used here. 
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Fig. 15. — Scaling exponents of the structure functions of column density obtained from 
densities shaped in the form of a cloud. Plots show exponents of the structure functions of 
first to sixth order excluding the third one, which is taken as a normalization. For comparison 
the dotted lines show the corresponding values of exponents for the original density, i.e. not 
shaped in the form of a cloud (the same as in Figure [8]). The low-resolution models are used 
here as well. 
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Fig. 16. — PDFs for density (left) and power spectra of density (right) for experiments with 
Ms ~ 0.7 and Ma ~ 0.7 for three different resolutions: 128^, 256^, and 512^. Spectra in 
the right plot have been shifted to make the dissipative parts coincide. 
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Fig. 17. — Scaling exponents for density for models with M.s ~ 0.7 and M.a ~ 0.7 for 
three different resolutions: 128'\ 256^, and 512^. Lines show theoretical scalings: dotted 
for Kolmogorov and dash-dotted and dashed for She-Leveque with parameter C = 1 and 2, 
respectively (fif = 3 and x = |). 
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Fig. 18. — Topology of the density structures (left column)^ and topology of the density 
integrated along the magnetic field {right column). Subsonic models have symmetric dis- 
tributions of the genus line, and the supersonic models have much longer tails of higher 
densities (see text for explanation). In all plots, vertical dashed lines designate the position 
of the mean density value. 
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Fig. 19. — Topology of the logarithm of density {left column) and the logarithm of column 
density integrated along the line of sight [right column). Plots correspond to models from 
Fig. [m In all plots, vertical dashed lines designate the position of the mean density value. 
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Table 1. List of All Analyzed Models 



Model Bext Pini -Ms -Ma Resolutioii Max. Time Number of Snapshots max Sp max 6V 



BIPIO 


1, 


.0 


10.0 


0.20±oo2 


0.62±oo4 


128^ 


B1P9 


1, 


.0 


9.0 


0.2i±oo2 


0j2±oo4 


128^ 


B1P8 


1, 


.0 


8.0 


0.22±oo2 


0.62±o-04 


128^ 


B1P7 


1, 


.0 


7.0 


0.24±oo2 


0.63±o°4 


128^ 


B1P6 


1, 


.0 


6.0 




0.63±o°^ 


128^ 


B1P5 


1, 


.0 


5.0 


0.28±°-°2 


0.63±""'* 


128^ 


B1P4 


1 


.0 


4.0 


32±o.o3 


64±"04 


128^ 


B1P3 


1, 


.0 


3.0 


0_37±o.o4 


0_64±oo5 


128^ 


B1P2 


1, 


.0 


2.0 


0.46±o.05 


0.65±oo5 


128^ 


BlPl 


1 


.0 


1.0 


oge±o.o5 


0gg±0.05 


128^ 


B1P.9 


1 


.0 


0.9 


0J0±O.O5 




128^ 


B1P.8 


1, 


.0 


0.8 


QJ4±0.06 




128^ 


B1P.7 


1, 


.0 


0.7 






128^ 


B1P.6 


1 


.0 


0.6 


0.86±o"^ 




128-3 


B1P.5 


1, 


.0 


0.5 


0.94±o.08 


0_66±oo5 


1283 


B1P.4 


1, 


.0 


0.4 


1 Q5±0.08 




1283 


B1P.3 


1, 


.0 


0.3 


1.21±0"9 


Q gg±0.05 


1283 


B1P.2 


1 


.0 


0.2 


;^ 4y±0.09 


0.66±°°^ 


1283 


BlP.l 


1, 


.0 


0.1 


2.1±0-2 


0.67±0-06 


1283 


B.IPIO 


0, 


.1 


10.0 


0.23±0-0i 


7 i±0.3 


1283 


B.1P9 


0, 


.1 


9.0 


24±" "i 


7 -[^±0.4 


1283 


B.1P8 


0, 


.1 


8.0 


25±o.02 


7 -|^±0.4 


1283 


B.1P7 


0, 


,1 


7.0 


Q_27±o.o2 


7 ]^±0.4 


1283 


B.1P6 


0, 


.1 


6.0 


0.29±oo2 


7_2±o.3 


1283 


B.1P5 


0, 


.1 


5.0 


32±o.o2 


7_2±o.3 


1283 


B.1P4 


0, 


.1 


4.0 


0.36±°02 


7_2±o.3 


1283 


B.1P3 


0, 


.1 


3.0 


0.42±o.02 


73±o.3 


1283 


B.1P2 


0, 


.1 


2.0 


0.52±oo3 


7.3±o-4 


1283 


B.lPl 


0, 


.1 


1.0 


0j2±o.03 


7_2±o.3 


1283 


B.1P.9 


0, 


.1 


0.9 


0.77±o.03 


7.3±0-3 


1283 


B.1P.8 


0, 


.1 


0.8 


0.82±0-03 


7_3±o.2 


1283 


B.1P.7 


0, 


.1 


0.7 


0,88±o.o3 


7 4±o.3 


1283 


B.1P.6 


0, 


.1 


0.6 


0.94±o.o3 


7.3±o-3 


1283 


B.1P.5 


0, 


,1 


0.5 


l_03±oo3 


7.3±0-3 


1283 


B.1P.4 


0, 


.1 


0.4 


1 15±0.04 


73±o.3 


1283 


B.1P.3 


0, 


.1 


0.3 


1 33±o.o4 


7.3±"-3 


1283 


B.1P.2 


0, 


.1 


0.2 


1.66±°06 


7_4±o.3 


1283 



25.0 


16 


~0.06 


—0.6 


25.0 


16 


~0.06 


—0.6 


25.0 


16 


~0.06 


—0.6 


25.0 


16 


~0.06 


—0.6 


25.0 


16 


^0.06 


—0.6 


25.0 


16 


^0.07 


—0.6 


25.0 


16 


~0.08 


—0.6 


25.0 


16 


~0.1 


—0.6 


25.0 


16 


~0.1 


—0.6 


25.0 


16 


~0.2 


—0.6 


25.0 


16 


~0.2 


—0.6 


25.0 


16 


~0.3 


—0.6 


25.0 


16 


~0.3 


—0.6 


25.0 


16 


~0.3 


—0.6 


25.0 


16 


~0.3 


—0.6 


25.0 


16 


^0.4 


—0.7 


25.0 


16 


~0.4 


—0.7 




16 

-LU 


r^O 5 


— n 7 


25.0 


16 


~0.9 


—0.7 


25.0 


16 


~0.06 


—0.7 


25.0 


16 


~0.06 


—0.7 


25.0 


16 


~0.06 


—0.7 


zo.U 


it) 


~U.Ut) 


— U. 1 


25.0 


16 


~0.06 


-0.7 


25.0 


16 


~0.06 


-0.7 


25.0 


16 


~0.07 


-0.7 


25.0 


16 


~0.09 


-0.7 


25.0 


16 


~0.1 


-0.6 


25.0 


16 


-0.2 


-0.7 


25.0 


16 


~0.3 


-0.7 


25.0 


16 


~0.3 


-0.7 


25.0 


16 


~0.3 


-0.7 


25.0 


16 


-0.3 


-0.7 


25.0 


16 


-0.4 


-0.7 


25.0 


16 


~0.4 


-0.7 


25.0 


16 


-0.5 


-0.7 


25.0 


16 


-0.6 


-0.8 
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Table 1 — Continued 



Model 


Bext 


ini 


Ms 


Ma 


Resolution 


Max. Time 


Number of Snapshots 


max 6p 


max 6V 


B.lP.l 


0.1 


0.1 


2.3±0-i 


7 4±o.3 


128^ 


25.0 


16 


~0.8 


-0.8 


B1P8 


1.0 


8.0 


0.23±""i 


0.64±oo2 


256^ 


10.0 


8 


~0.05 


-0.6 


B1P4 


1.0 


4.0 


0.33±"oi 


0.65±oo2 


2563 


15.5 


14 


-0.1 


-0.7 


BlPl 


1.0 


1.0 




0_68±oo3 


2563 


10.1 


6 


-0.3 


-0.6 


BlP.l 


1.0 


0.1 


2.20±oo9 


0.69±oo3 


2563 


10.1 


9 


~0.9 


-0.7 


BlP.Ol 


1.0 


0.01 




OJO±0.04 


2563 


10.1 


10 


-1.9 


-0.7 


B.1P8 


0.1 


8.0 


0.26±oo3 


7.3±o-7 


2563 


14.5 


14 


-0.03 


-0.7 


B.1P4 


0.1 


4.0 


0.36±oo4 


7_2±o.7 


2563 


16.7 


14 


-0.07 


-0.7 


B.lPl 


0.1 


1.0 


Q y4±0.06 


7_4±o.6 


2563 


21.0 


16 


-0.2 


-0.7 


B.lP.l 


0.1 


0.1 


2 34±o.o8 


7 4±o.3 


2563 


22.3 


17 


-0.8 


-0.7 


B.lP.Ol 


0.1 


0.01 




71±0.3 


2563 


20.9 


16 


-1.5 


-0.7 


BlPl 


1.0 


1.0 


-0.6 


~0.6 


5123 


7.2 


1 


-0.3 


-0.6 
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Table 2. Slopes of the Power Spectrum of Density and the Logarithm of Density 

Fluctuations 





Ma ~ 0.7 










Ms 


ftp 




Ms 


ftp 


'-'^logp 


0.23±o.oi 


-2.3±o-3 


-2.3±o-3 


0,26±o.o3 


_1 7±0.3 


_1 7±0.3 


0.33±°-oi 


_2_2±o.3 


_2_2±o.3 


0.36±o-o^ 


.1 7±0.3 


.1 7±0.3 


QQ8±om 


.2.o±o.3 


_2.i±o.3 


Q 74±o.o6 


_IQ±0.2 


.1 Q±0.3 


2 20±o-03 


-1.3±o-2 


.2.o±o.2 


2_34±o.o8 


-1.2±o-2 


_l_g±0.2 


70±o.3 


-0.5±o-i 


.l_7±o.2 


7 -|^±0.3 


-0.6±°-2 


-1.5±°-2 



Note. — The values have been estimated within the inertial 
range for models with Ma ~ 0.7 (left) and Ma ~ 7 (right). 
Errors of spectral indices combine the errors of estimation at 
each time snapshot and the standard deviation of variance in 
time. Errors for sonic Mach numbers are the standard deviation 
of their variance in time calculated over the period starting from 
t > 5 to the last available snapshot. 



